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The stability of a Bose-Einstein condensed state of trapped ultra-cold atoms is investigated under 
the assumption of an attractive two-body and a repulsive three-body interaction. The Ginzburg- 
Pitaevskii-Gross (GPG) nonlinear Schrodinger equation is extended to include an effective potential 
dependent on the square of the density and solved numerically for the s— wave. The lowest fre- 
quency of the collective mode is determined and its dependences on the number of atoms and on 
the strength of the three-body force are studied. We show that the addition of three-body dynamics 
can allow the number of condensed atoms to increase considerably, even when the strength of the 
three-body force is very small compared with the strength of the two-body force. We also observe a 
first-order liquid-gas phase transition for the condensed state up to a critical strength of the effective 
three-body force. 
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I. INTRODUCTION 

The theoretical research on Bose-Einstein condensa- 
tion (BEG) a phenomenon predicted more than 70 
years ago, is receiving considerable experimental and the- 
oretical support in recent years j^. The relevance of 
BEG for understanding the properties of liquid ^He vi^as 
pointed out by London suggesting that the peculiar 
phase transition that liquid helium undergoes at 2.18K is 
a BEG phenomenon. It is also important to observe that, 
at the level of two-body collisions, Bogoliubov in 1947 Q 
has shown for homogeneous gas that BEG is only possible 
for systems with repulsive potentials. 

Intense experimental researches on BEG for magneti- 
cally trapped weakly interacting atoms have been done 
recently [^-||. In the experiment reported in 0, a 
condensate of approximately 2000 spin-polarized °^Rb 
atoms was produced in a cylindrically symmetric mag- 
netic trap (2|j9[]. It is a common understanding that, 
in low temperature and density, where interatomic dis- 
tances are much greater than the distance scale of atom- 
atom interactions, two-body interactions take a simple 
form, and three-body interactions can be neglected. At 
such regime, only two-body s— wave scattering is impor- 
tant. With temperature low enough the magnitude of the 
scattering length a is much less than the corresponding 
thermal de Broglie wavelength, and the exact shape of 
two-atom interaction is unimportant. 

The experimental evidences of Bose-Einstein conden- 
sation (BEG) in magnetically trapped weakly interacting 
atoms ||-^ brought a considerable support to the theo- 
retical research on bosonic condensation. The nature of 
the effective atom-atom interaction determines the stabil- 
ity of the condensed state: the two-body pseudopoten- 



tial is repulsive for a positive s— wave atom-atom scat- 
tering length and it is attractive for a negative scattering 
length . The ultra-cold trapped atoms with repulsive 
two-body interaction undergoes a phase-transition to a 
stable Bose condensed state, in several cases found exper- 
imentally, as for s^Rb @, ^SNa and |]. However, 
a condensed state of atoms with negative s— wave atom- 
atom scattering length (as in case of ''Li Q) would be 
unstable, unless the number of atoms N is small enough 
such that the stabilizing force provided by the harmonic 
confinement in the trap overcomes the attractive inter- 
action, as found on theoretical grounds ||ll|,|l^. It was 
indeed observed in the ^Li gas for which the s— wave 
scattering length is a = —14.5 ± 0.4 A, that the number 
of allowed atoms in the Bose condensed state was lim- 
ited to a maximum value between 650 and 1300, which 
is consistent with the mean- field prediction pl] |. 

So, for systems of atoms with attractive two-body in- 
teraction, it is widely believed |11 l^Jl^ that the conden- 
sate has no stable solution above certain critical number 



of atoms N„ 



However, in this case the addition of 



a repulsive potential derived from three-body interac- 
tion is consistent with a number of atoms larger than 
Nmax- Even for a very small strength of the three-body 
force, the region of stability for the condensate can be 
extended considerably, as previously reported in p^ , p^ , 
and shown in more detail in the present work. By consid- 
ering the possible effective interactions, it was reported 
in Ref. ||l^ that a sufficiently dilute and cold Bose gas 
exhibits similar three-body dynamics for both signs of 
the s— wave atom-atom scattering length. It was also 
suggested that, for a large number of bosons the three- 
body repulsion can overcome the two-body attraction, 
and a stable condensate will appear in the trap |jl^. If 
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an atomic system is characterized by having efTectively 
an attractive two-body interaction together with a repul- 
sive three-body interaction, two mechanisms for stabihty 
are possible: (a) the kinetic energy acting at lower densi- 
ties and (b) the repulsive weak three-body force effective 
at higher densities. These mechanisms indicate that, for 
the same number of atoms, one lower-density phase and 
a higher-density phase can be found, if the three-body 
force is weak enough not to dominate the effective inter- 
action. 

It was pointed out in Rcf. p9[ | that an easier experi- 
mental approach to probe density fluctuations is to con- 
sider an observable directly sensitive to the probability 
of finding three atoms near each other, which will cor- 
respond to the loss rate of atoms due to three-body re- 
combination. Such a three-body recombination rate in 
BEC, was considered recently in Refs. jH], ||2l| and |14] 
(see also the review of [|3|). It was shown in Ref. |2C] 
that the three-body recombination coefficient of ultra- 
cold atoms to a weakly bound s level goes to infinity 
in the Efimov limit p^ . The Efimov limit is a particu- 
larly interesting three-body effect, which happens when 
the two-body scattering length is very large (positive or 
negative). In this case, with the two boson energy close 
to zero, the three-boson system presents an increasing 
number of loosely bound three body states, which have 
large spatial extension and do not depend on the details 
of the interaction |23|. So, our main motivation here is to 
provide an extension to the GPG equation |^,^, which 
considers a three-body interaction and, in this way, pro- 
vides the framework for a numerical investigation of the 
relevance of three-body interaction in Bose-Einstein con- 
densation. 

In the present work we consider a possible general sce- 
nario of atomic systems with attractive two-body and re- 
pulsive three-body interactions. By using the mean-field 
approximation, we investigate the competition between 
the leading term of an attractive two-body interaction, 
originated from a negative two-atom s— wave scattering 
length, and a repulsive three-body effective interaction, 
which can happen in the Efimov limit (|a| oo) as 
discussed in Ref. jl^] Q We show that, in a dilute gas, 
a small repulsive three-body force added to an attractive 
two-body interaction is able to stabilize the condensate 
beyond the critical number of atoms in the trap, found 



just with attractive two-body force |11|, such that a kind 
of liquid-gas phase-transition occurs. The plan of the pa- 
per is as follows. In section II, we introduce the Ginzburg 
- Pitaevskii - Gross (GPG) formalism. In section III, we 
present the main numerical results for the static solu- 
tions, together with a variational analysis. In section IV, 



*The physics of three-atoms in the Efimov limit is discussed 
in Ref. |2d] , that extends a previous study of universal aspects 
of the Efimov effect [El. 



we present a stability analysis and results for collective 
excitation in the condensate. In this section IV we also 
observe that the inclusion of three body effects points out 
possible evidences of a liquid - gas phase transition in the 
condensate. Finally, in section V, we present our main 
conclusions. 



II. GINZBURG - PITAEVSKII - GROSS 
FORMALISM 

In the following, we present our formalism, where the 
original Ginzburg - Pitaevskii - Gross (GPG) non-linear 
equation ||2^,|2^, which includes a term proportional to 
the density (two-body interaction), is extended through 
the addition of a term proportional to the squared- 
density (three-body interaction). Next, after reducing 
such equation to dimensionless units, we study numeri- 
cally the s— wave solution by varying the corresponding 
dimensionless parameters, which are related to the two- 
body scattering length, the strength of the three-body 
interaction and the number of atoms in the condensed 
state. As particularly observed in Ref. |2^, to incorpo- 
rate all two-body scattering processes in such many par- 
ticle system, the two-body potential should be replaced 
by the many-body T— matrix. Usually, at very low en- 
ergies, this is approximated by the two-body scattering 
matrix, which is directly proportional to the scattering 
length So, in order to obtain the desired equation, 
we first consider the effective Lagrangian, which describes 
the condensed wave-function in the Hartree approxima- 
tion, implying the GPG energy functional: 
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In our description, the atomic trap is given by a rota- 
tionally symmetric harmonic potential, with angular fre- 
quency Lo, and £i gives the effective atom interactions up 
to three particles. 

The effective interaction Lagrangian for ultra-low tem- 
perature bosonic atoms, including two- and three-body 
scattering at zero energy, is written as: 

A = J dVidV2dVidV^«'t(Pi)*t(p2)*(ri)*(r2) 
X (r^i2 |r(2)(0)| fi2) (rl + r'2 - r i - rl) 
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where r'12 and ^3 are the relative coordinates, given 
by fi2 = ri — f2 and ^3 = "^3 — (^1 + ^2)/2; and 
^123 = (fi + fa + f3) . T(3) (0) and ^ (0) are the corre- 
sponding three-body T— matrix and two-body T— matrix 
for the pair jk, which are evaluated at zero-energy. The 
two-body T— matrix for each pair (jk) is subtracted from 
T^^^ (0) to avoid double counting and Ki is the kinetic en- 
ergy operator for particle i. 

We can approximate the above effective interaction La- 
grangian at low densities by averaging the T— matrices 
over the relative coordinates, considering that the ther- 
mal wave-length is much greater than the characteristic 
interaction distances. 
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The integrations of the T-matrices over the relative co- 
ordinates gives the zero momentum matrix elements: 



d^r[2d^n2(r'i2 



r(2)(0) 



f^l2^ = 



(27r)3/pi2-0 r(2)(0) 



P12 = 



(4) 



where a is the two-body scattering length. For the con- 
nected three-body T— matrix, also by integrating over 
the coordinates, we obtain the corresponding zero mo- 
mentum (pi2 = 0, P3 = 0) matrix elements, which give 
us the strength of the three-body effective interaction A3, 
as follows: 



J d^r[2d^Ri^d^ri2d^R3 
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= 2X. 
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where (| = (pi2 = 0, P3 = 0| and |) = \pi2 = 0, F3 = 0). 

The nonlinear Schrodinger equation, which describes 
the condensed wave-function in the mean-field approxi- 
mation, is obtained from the effective Lagrangian given 
in Eq. (|^). By considering the interaction in Eq. (g), it 
can be written as [ESl 
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+ A3iV2|vI/(f,i)|4]vI/(f,i). 
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For a stationary solution, ^{f, t) = e ^t^*/'^ '^(r), and the 
above equation can be written as 



2m 2 



m 



(7) 

where fi is the chemical potential (single particle energy) 
and ip{r) is normalized as 

dV|V'(f)|^ = 1. (8) 
The total energy of the system is given by 

J 2m 2 

- ' + ^A3hA r ' . 9 

2 m 3 J 

The central density of the system can be obtained di- 
rectly from the solution of the above equation, normal- 
ized as in Eq. (||): 

Pc = iV|V'(0)p. (10) 

The physical scales presented in the above equations can 
be easily recognized by working with dimensionless equa- 
tions. By rescaling Eq. (^ for the s— wave solution, we 
obtain 



d' 1 2 
dx^ 4 



93- 



$(a;) = /3$(x) , 
(11) 



where x = ^J2mwjfi r and $(x) = N^^"^ ^Ji■n\a\ rip{r). 
The dimensionless parameters, related to the chemical 
potential and the three-body strength are, respectively, 
given by 

(12) 



P = - — and 33 = A3?ia; 



The normalization for <i>(x), obtained from Eq. (g), de- 
fines a number n related to the number of atoms N: 



2mu} 



/ dx\^{x)\'^ = n, where n = 2N\a\ 
Jo 

The boundary conditions in Eq.(|l^) are given by |^l[ 
$(a; ^ 0) 



(13) 



^{x —> 00) cx exp — 



ln(a;) 



(14) 



In terms of the dimensionless variables, the total energy 
of the system is given by 



E = huiN I dx ■ 



d(j){x) 



dx 



X (j) (x) 



n(j)'^{x) n^g3(t)^{x) 



2a;2 3.t4 
where (f>{x) = $(x)/n^/^ is normalized to one. 



(15) 
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III. LIQUID-GAS PHASE TRANSITION - STATIC 
SOLUTIONS 

A. Variational Approach 

As a further reference to our results, and the stabihty 
analysis, it will be helpful first to consider a variational 
procedure using a trial gaussian wave function for 
tp{'r). So, in Eq. we consider the following variational 
wave function (normalized to one): 



a 



Ipvarir) = 



1 muj\ * 
— — exp 



r f muj \ 
2^ VY) 



(16) 



where a is a dimensionless variational parameter. The 
corresponding root-mean-square radius, tq, will be pro- 
portional to the variational parameter a: 



ro = V {r'^)var = a 



3h 

2muj 



(17) 



The expression for the total variational energy, which is 
obtained after replacing Eq. (|l^) in Eq. (||), is given by 



9\/37i 



(18) 



In the same way, we can obtain the corresponding varia- 
tional expression for the single particle energy, Eq. (|^): 



(19) 



The variational central density, using Eqs. (|l^) and ^ 
can also be given in terms of this parameter a: 



ttTi / a'^ 



(20) 



The approximate solutions for the total energy are ob- 
tained from the extrema of ([l^) with respect to variation 
of the parameter a. 

The variational solutions of E^ar{ct) are given, as a 
function of n and (where a < and 93 > 0), by the 
real roots of dEyar[o)/ da = |. 



^By using a numerical procedure one can reach easily the 
extrema of Eq. ( |l^ by varying the parameter a, once the 
other parameters are fixed. 




FIG. 1. In the lower part, we have a comparison between 
variational (solid curve) and exact (dashed curve) numerical 
calculations of the condensate energy as a function of the 
reduced number of atoms n for 33 = 0.005 . In the upper 
frame we show five plots of the variational energy as a function 
of the variational parameter a for five particular values of 
n shown also in the lower frame. (I) (resp IV) corresponds 
to a small (large) n region where only one stable solution is 
encountered; (II) (resp III) to a small (large) n region where 
we observe three extrema for the energy; (C) corresponds to 
a particular n where we obtain two stable solutions with the 
same energy E\ = E2- E is given in units of {Nhu))/n. 

In Fig. 1, we first illustrate the variational procedure 
considering an arbitrarily small three-body interaction, 
chosen as 33 = 0.005. In the upper part of the figure, 
we show five small plots for the total variational energy 
E, in terms of the variational width a. Each one of the 
small plots corresponds to particular values of n. For 
each number n we report the energy of the variational 
extrema in the lower part of Fig. 1. In region (I) where 
the number of atoms is still small, the attractive two 
body force dominates over the repulsive three-body force 
and just one minima of the energy as a function of the 
variational parameter a is found. That is also the case 
for 33 — 0. When the number of atoms is further in- 
creased (region (II)) two minima appear in the energy 
E (a) . An unstable maximum is also found between the 
two minima. The lower energy minimum is stable while 
the solution corresponding to the smaller a is metastable. 
This solution has a higher density and, consequently, its 
metastability is justified by the repulsive three-body force 
acting at higher densities. The minimum number n for 
the appearance of the metastable state is characterized 
by an inflection point in the energy as a function of a. 
The value of n at the inflection point corresponds to the 
beak in the plot of extremum energy versus n because 
for larger n three variational solutions are found as de- 
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picted in the lower part of Fig. 1. The attractive two- 
body and trap potentials dominate the condensed state 
in the low-density stable phase up to the crossing point 
(C). At this point, the denser metastable solution be- 
comes degenerate in energy with the lower-density stable 
solution and a first order phase transition takes place. 
Since the two solutions differ by their density this tran- 
sition is analogous to a gas-liquid phase transition for 
which the density difference between the liquid and the 
gas is the order parameter. In the variational calcula- 
tion this occurs at the transition number n «1.3 while 
the numerical solution of the NLSE gives 1.2. In region 
(III), we observe two local minima with different ener- 
gies, a higher-density stable point and a lower-density 
metastable point. The metastable solution disappears in 
the beak at the boundary between region (III) and (IV). 
In regions (III) and (IV) the three-body repulsion sta- 
bilized a dense solution against the collapse induced by 
the two-body attraction. The qualitative features of the 
variational solution is clearly verified by the numerical 
solution of the NLSE, as shown by the dashed curve. 



B. Numerical Results 



The numerical solutions of Eq. (y_l|) are obtained for 
several values of /3, using three values of 173 to char- 
acterize the solutions. We have used the Runge-Kutta 
(RK) and "shooting" method to obtain the correspond- 
ing solutions in each case po|| . The stability assignment 
for the stationary solutions was made by studying the 
corresponding time dependent Schrodinger equation, us- 
ing the Crank-Nicolson (CN) method (see Refs. |ll[] and 
|31| ) . The numerical procedure to determine such stabil- 
ity was done in the following way: when applying the CN 
method, we started by using the static solution obtained 
from the RK method and observed if the modulus of the 
wave function remained constant. If this was occurring 
for a long period of time (of about 500 units of dimen- 
sionless time r = cot) the solution was considered stable, 
otherwise unstable. 

In Fig. 2 we present the total energy as a function of 
the number of atoms, represented by the reduced number 
n defined in Eq. (13), for three significative values of the 
quintic parameter (73, given by 0, 0.016 and 0.03. The 
results agree with Ref. [|l5|. When (73 = 0, the stable 
solutions for the energy starts at zero (for n = 0) and 
reaches a critical limit at Umax — 1-62. There is no solu- 
tions for higher n, but the plot also shows a branch with 
unstable solutions (with higher energies) for n < 1.62. 
Our results are consistent with results given in Ref. . 
When (73 = 0.03, only stable solutions appear for the 
energy, with no limit in the number of atoms, having a 
maximum at n 2. So, this and higher values for (73 
already represent a dominance of the quintic term in the 
interaction of Eq. (11). We observe that the numerical 
stability analysis is consistent with the variational ap- 



proach discussed in the previous sub-section. The more 
interesting case represented in Fig. 2 is for (73 =0.016, 
as in such a case we observe a region of the plot where 
we can have up to three solutions for the same n. The 
inset to this figure amplifies the region of the plot where, 
for 33 — 0.016, the solutions become unstable (between 
A and B) or metastable (between A and C, or B and 
C). At the point C a phase transition occurs from a less 
denser (gas) to a more denser (liquid) phase. 



93=0 

93=0.016 

93=0.03 
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n 

FIG. 2. The total energy, in 

units of {Tiuj)/ is shown as a function of 

the reduced number of atoms n, given by Eq. (^3|), for 33 =0, 
0.016 and 0.03. The inset points out critical limits discussed 
in the text. 

In Fig. 3, following a correspondence to Fig. 2, we 
present the results for the chemical potential in dimen- 
sionless units (/3) as a function of n. The line with arrow 
in the inset to this figure indicates the approximate posi- 
tion in n, where the phase-transition (from a 'gas' phase 
to a 'liquid' phase) occurs. For 53 = 0.016 the part of the 
plot linking points A and B is unstable (see both Figs. 
2 and 3), otherwise it is stable. Finally, for = 0.03, 
the function of the energy in terms of n is always single 
valued and stable. Our calculation for 93 = also agrees 
with results presented in Ref. [lll| , with the maximum 
number of atoms limited to rimax ~ 1.62 0. As we can 
see, for n < Umax two solutions are possible, one of them 
being unstable. For 173 higher than zero, a new pattern 
appears. For instance, the plot for the case of 33 = 0.016 
(see the inset) can be divided in several sectors according 
to the stability analysis, with the help of Fig. 2: Start- 
ing from 71 = (/3 = 1.5) until point Cq, and from Cl 



tQur n is equal to \Clf\ of Ref. (Tl). 
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to higher values of n, we have stable solutions; from Cg 
to A and from B to Cl we obtain metastable solutions; 
from A to B the solutions are unstable, corresponding to 
maxima for the energies. 
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g, = 0.016 
g, = 0.03 
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FIG. 3. The chemical potential, in dimensionless units 
{13 = fi/{hu!), is shown as a function of the reduced num- 
ber of atoms n, for the same set of gs shown in Fig. 2. The 
inset points out the critical limits corresponding to Fig. 2 (Cg 
and Cl corresponds to C), and the straight line with arrow 
indicates the transition from a less denser to a more denser 
phase. 

In Fig^ we also plot the central density pc, defined 
in Eq. (10), as a function of the number n. We use the 
same values of the parameter as used in Figs. 2 and 
3. The inset to the figure also points out the phase tran- 
sition which occurs when =0.016. As the straight 
line with arrow shows, after the transition the system 
becomes more than three times denser than the original 
one. Also, for < 33 < 0.0183, we observe that the 
density pc presents back bending typical of a first order 
phase transition. 

By extending the observations of a first order phase 
transition, given in Figs. 2-4 for 173 =0.016, we also 
determined the region of gs where such kind of phase- 
transition can occur. In Fig. 5 we have a phase-diagram, 
where it was shown the critical boundary separating 
the two phases and a critical point at n = 1.8 and 
(73 = 0.0183. For g^ less then such critical value, we 
observe two regions with distinct phases, similar to gas 
and liquid phases. These two different phases are also 
clearly identified in our Fig. 4, where we present the 
central density as a function of n. 
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FIG. 4. Central density, in dimensionless units, as a func- 
tion of the number n, for the same set of parameters gs given 
in Figs. 2 and 3. 
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FIG. 5. Graphical representation of the interface of the two 
distinct phase (gas and liquid), in the plane defined by the re- 
duced number of atoms n and the parameter gs (lower frame); 
and for the central density pc versus pa (upper frame). The 
arrows in the lower frame correspond to the point where it 
occurs the phase transition for gg =0.016, when changing n. 

For each g^ , the transition point given by the crossing 
point in the E versus n (see Fig. 2) corresponds to a 
Maxwell construction in the diagram of p versus n. At 
this point an equilibrated condensate should undergo a 
phase transition from the branch extending to small n 
to the branch extending to large n. The system should 
never explore the back bending part of the diagram be- 
cause, as seen in Fig. 2, it is an unstable extremum of the 
energy. From Figs. 1-5, it is clear that the first branch 
is associated with small densities, large radii, and posi- 
tive chemical potentials while the second branch presents 
a more compact configuration with a smaller radius a 
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larger density and a negative chemical potential. This 
justify the term gas (G) for the first one and liquid (L) 
for the second one. However we want to stress that both 
solutions are quantum fluids. 



IV. COLLECTIVE EXCITATIONS 

In this section, from the time evolution of GPG equa- 
tion, given in Eq. (^), we consider the ground-state col- 
lective excitations for the system [^2|-p4[. Following 
Ref. iQ, the collective excitations are described by the 
Bogoliubov equations [H,p5 29 35|] . After including three 



body interactions they take the form 



+ hLO,]v, + {NUo + 2\sN^\i;g\^}[rfu, = 0, (21) 



{NUo + 2\3N'\ijg\'}[ijgYv, 
^.1 ..... 



where 



C,=Ho-^i + 2UoN + 3X3N^\^g\\ (22) 

Hq is the harmonic oscillator hamiltonian, Uq = 
— (47r?i^|a|)/m, is the frequency of the collective oscil- 
lations, N is the number of atoms and tpg = ipg{v) is the 
ground state solution of the Eq. (|^) . The above equations 
have been solved by using several methods [^^4|j3^ . In 
the present calculations we have employed two methods: 
a time-dependent and a time-independent one. In the 
time-dependent procedure we have added a weak pertur- 
bation to the potential and, with CN algorithm, exam- 
ined the time evolution of Eq. (^) for a selected point of 
the wave-function. The lowest collective oscillations (wi^) 
were determined through Fourier transformation . By 
using the time independent algorithm, we have solved 
Eqs. (|2^) with the matching algorithm |37j generalized 
for two functions u and v. The method works by de- 
parting from the analytically known u, v and for the 
harmonic oscillator (chemical potential near to 3/2tiLu). 
Then we successively apply the matching method for the 
coupled u and v, gradually decreasing the chemical po- 
tential. This allows to reach subsequent solutions, by em- 
ploying the deformation algorithm described in Ref. . 
We obtain exact agreement between both methods, time- 
dependent or time-independent one. 

Figure 6 shows the collective frequencies cui, as a func- 
tion of n for the first mode {1 = 0). The solutions cor- 
responding to 93 = agree well with the ones given in 
Ref. loosing stabihty as lUi, ^ 0. By using this 

critcrium, we have obtained the regions of stability for 
53 = 0.016. For 33 = 0.03 all the solutions are stable. 
Following the inset of Fig. 6, for = 0.016, one can 
observe that, as the number of atoms is increased, in the 
less denser phase, the frequency of the collective exci- 
tations decreases and are related to stable solutions till 
the point Cg; from this point down to the point A (in- 
creasing n), the frequency continues to decrease to zero, 
but now related to meta-stable solutions. As already ex- 
plained previously in the discussion of Figs. 2-5, and also 



from the variational energy solutions given in Fig. 1, it is 
very likely that occurs a phase transition, from Cg to Cl 
(or from the meta-stable solutions, given in the branches 
Cg— A and B— Cl, to the corresponding stable solutions 
with fixed n). Once in the denser phase (from B pass- 
ing through the point Cl), the frequency of the collective 
excitations increases as the number of atoms increases, 
contrary to the behavior observed for the system in the 
less denser phase. This can be qualitatively understood 
considering the variational energy of the two phases and 
the corresponding stable energy as shown in Fig. 1. The 
curvature of the variational energy as a function of a at 
the minimum for the liquid phase is bigger than the cor- 
responding one in the gas phase [compare in Fig.l the 
insets (I) and (II) with the insets (III) and (IV)]. This 
indicates, in agreement with Fig. 6, that the restoration 
force is stronger for the liquid phase than for the gas 
phase and consequently the frequencies of the collective 
modes starting at the point Cl are higher than the cor- 
responding ones for the gas phase ending at Cg- As we 
include more particles the frequencies of the oscillations 
increase in the liquid phase. Corresponding to Fig. 6, 
in Fig. 7 the collective frequencies are shown as a func- 
tion of the chemical potential /3. From right to left, as 
the chemical potential decreases till Cg, (3 also decreases; 
from Cg to A, and from B to Cl the solutions are meta- 
stable, such that the system will look for a transition to 
a stable branch (from Cg, increasing /3, and from Cl, de- 
creasing /3). From Cl, as we further decreases the value 
of P the frequency of the collective excitations increases. 
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FIG. 6. Collective frequencies as a function of the reduced 
number of atoms n. The inset shows the critical points cor- 
responding to the previous figures. 
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FIG. 7. Collective frequencies as a function of the chemical 
potential /9. 

V. CONCLUSIONS 

To summarize, we have presented results for the total 
energy, chemical potential, central density, in terms of 
the number of atoms in the condensed state, for a range 
of values of the three-body strength. We also study the 
lowest collective mode excitations of the ground-state. 
Our calculation presents, at the mean-field level, the con- 
sequences of a repiilsive three-body effective interaction 
for the Bose condensed wave-function, together with an 
attractive two-body interaction. A first-order liquid-gas 
phase-transition is observed for the condensed state as 
soon as a small repulsive effective three-body force is 
introduced. In dimensionless units the critical point is 
obtained when g-^ « 0.0183 and n « 1.8. The characteri- 
zation of the two-phases through their energies, chemical 
potentials, central densities and radius were also given 
for several values of the three-body parameter ^3. The 
results presented in this paper can be relevant to deter- 
mine a possible clear signature of the presence of repul- 
sive three-body interactions in Bose condensed atoms. It 
points to a new type of phase transition between two 
Bose fluids. Because of the condensation of the atoms in 
a single wave-function this transition may present very 
peculiar fluctuations and correlations properties. As a 
consequence, it may fall into a different universality class 
than the standard liquid-gas phase transition, which are 
strongly affected by many-body correlations. This mat- 
ter certainly deserves further studies. 
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